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ABSTRACT 

We describe local shearing box simulations of turbulence driven by the magnetorotational instability (MRI) 
in a collisionless plasma. Collisionless effects may be important in radiatively inefficient accretion flows, 
such as near the black hole in the Galactic Center. The MHD version of ZEUS is modified to evolve an 
anisotropic pressure tensor A fluid closure approximation is used to calculate heat conduction along magnetic 
field lines. The anisotropic pressure tensor provides a qualitatively new mechanism for transporting angular 
momentum in accretion flows (in addition to the Maxwell and Reynolds stresses). We estimate limits on the 
pressure anisotropy due to pitch angle scattering by kinetic instabilities. Such instabilities provide an effective 
"collision" rate in a collisionless plasma and lead to more MHD-like dynamics. We find that the MRI leads to 
efficient growth of the magnetic field in a collisionless plasma, with saturation amplitudes comparable to those 
in MHD. In the saturated state, the anisotropic stress is comparable to the Maxwell stress, implying that the 
rate of angular momentum transport may be moderately enhanced in a collisionless plasma. 
Subject headings: accretion, accretion disks — MHD — methods: numerical — plasmas - turbulence 



1. INTRODUCTION 

Following the seminal work of Balbus & HawlevI il991l) . 
numerical simulations have demonstrated that magnetohy- 
drodynamic (MHD) turbulence initiated by the magnetorota- 
tional instability (MRI) is an efficient mechanism for trans- 
porting angular momentum in accretion disks (e.g., Hawley, 
Gammie, & Balbus 1995, hereafter HGB; see Balbus & Haw- 
ley 1998, for a review). For a broad class of astrophysical ac- 
cretion flows, however, the MHD assumption is not directly 
applicable. In particular, in radiatively inefficient accretion 
flow (RJAF) models for accretion onto compact objects, the 
accretion proceeds via a hot, low density, collisionless plasma 
with the pr oton temperature larger than t he electron temper- 
ature (see Nara van et. al.lll 9981 lOu ataert "2003. for reviews). 
In order to maintain such a two-temperature flow the plasma 
must be collisionless, and there are many cases where the 
Coulomb mean-free path is many orders of magnitude larger 
than the system size. Motivated by the application to RIAFs, 
this paper studies the nonlinear evolution of the MRI in a col- 
lisionless plasma, focusing on local simulations in the shear- 
ing box limit. 

Quataert, Dorland, & Hammett (2001; hereafter'ODHi) and 
Sharma, Hammett, & Quataert (2003; hereafter SHO) showed 
that the linear dynamics of the MRI in a collisionless plasma 
can be quite different from that in MHD. The maximum 
growth rate is a factor of ~ 1.7 larger and, perhaps more im- 
portantly, the fastest growing modes can shift to much longer 
wavelengths, giving direct amplification of long wavelength 
modes. Dynamical instability exists even when the magnetic 
tension forces are negligible because of the anisotropic pres- 
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sure response in a collisionless plasma. In related work using 
Braginskii's anisotropic viscosity, Balbus (2004) has termed 
this modification of the MRI in the lo w coUisionali ty limit the 
"magnetoviscous" instability (see also llslam & Balbus..2005i) . 

In this paper, we are interested in simulating the dynamics 
of a collisionless plasma on length-scales (~ disk height) and 
time-scales orbital period) that are very large compared to 
the microscopic plasma scales (such as the Larmor radius and 
the cyclotron period). Since the ratio of the size of the ac- 
cretion flow to the proton Larmor radius is ^ 10** for typical 
RIAF models, direct particle methods such as PIC (particle in 
a cell), which need to resolve both of these scales, are compu- 
tationally challenging and require simulating a reduced range 
of scales. Instead we use a fluid-based method to describe 
the large-scale dynamics of a collisionless plasma ("kinetic 
MHD," described in §2). The key differences with respect 
to MHD are that the pressure is a tensor rather than a scalar, 
anisotropic with respect to the direction of the local magnetic 
field, and that there are heat fluxes along magnetic field lines 
(related to Landau damping and wave-particle interactions). 
The drawback of our fluid-based method is, of course, that 
there is no exact expression for the heat fluxes if only a few 
fluid moments are retained in a weakly collisional plasma (the 
"closure problem"). We use results from Snyder, Hammett, & 
Dorland (1997; hereafter SHD) who have derived approxima- 
tions for the heat fluxes in terms of nonlocal parallel tempera- 
ture and magnetic field gradients. These heat flux expressions 
can be shown to be equivalent to multi-pole Pade approxi- 
mations to the Z-function involved in Landau damping. This 
approach can be shown to converge as more fluid moment s 
of the distribution function are kept ('Ham mett et alJfT993ft . 
just as an Eulerian kinetic algorithm converges as more grid 
points in velocity space are kept. These fluid-based meth- 
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ods have been appHed with reasonable success to modeling 
collisionless turbulence in fusion plasmas, generally com- 
ing within a factor of 2 of more complicated kinetic calcu- 
lations in st rong tu rbulence regimes (e.g., Dimits et al. 2000; 
|Pai-ker et alJ ll994t iHammett et al. 1993; Scott 2005), though 
there can be larger differences in weak turbulence regimes 
(see|Hammett et al. 1993; Dimits et al. 2000, and references 
therein). The simulations we report on here use an even sim- 
pler local approximation to the heat flux closures than those 
derived in SHD. While not exact, these closure approxima- 
tions allow one to begin to investigate kinetic effects with rel- 
atively fast modifications of fluid codes. 

In a collisionless plasma the magnetic moment /i is an adi- 
abatic invariant. Averaged over velocity space, this leads 
to conservation of (fi) = p±/(pB). As a result, pressure 
anisotropy with p± > /?|| is created as the MRI amplifies the 
magnetic field in the accretion flow. This pressure anisotropy 
creates an anisotropic stress (like a viscosity!) which can be 
as important for angular momentum transport as the magnetic 
stress (as we show below). The pressure anisotropy cannot, 
however, grow without bound because high frequency waves 
and kinetic microinstabilities feed on the free energy in the 
pressure anisotropy, effectively providing an enhanced rate 
of collisions that limit the pressure tensor anisotropy (lead- 
ing to more MHD-like dynamics in a collisionless plasma). 
We capture this physics by using a subgrid model to restrict 
the allowed amplitude of the pressure anisotropy. This sub- 
grid model (described in §2.3) is based on existing linear and 
nonlinear studies of instabilities driven by pressure anisotropy 
(e.g- lHasegawa 1969; Garv et al. 1997). 

The remainder of this paper is organized as follows. In the 
next section (§2) we describe Kulsrud's formulation of kinetic 
MHD (KMHD) and our closure model for the heat fluxes in 
a collisionless plasma. We also include a linear analysis of 
the MRI in the presence of a background pressure anisotropy 
and describe limits on the pressure anisotropy set by kinetic 
instabilities. In §3 we describe our modifications to the ZEUS 
code to model kinetic effects. §4 presents our primary re- 
sults on the nonlinear evolution of the MRI in a collisionless 
plasma. In §5 we discuss these results, their astrophysical im- 
phcations, and future work required to understand the global 
dynamics of collisionless accretion disks. 

2. GOVERNING EQUATIONS 

In the limit that all fluctuations of interest are at scales 
larger than the proton Larmor radius and have frequencies 
much smaller than the proton cyclotron frequency, a collision- 
less plasma can be described by the following magnetofluid 
equations (e.g., Kulsrud 1983; SHD) : 

dp 
(9V 



■ + V-(pV) = 0, 



(1) 



p— + p(V-V)V: 

as 

— = Vx(VxB), 
at 



(V X B) X B 
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(3) 
(4) 



P = p±I+ {p\i -pi_) bb = pi_I + U, 

where p is the mass density, V is the fluid velocity, B is the 
magnetic field, Fg is the gravitational force, b = B/ |B | is a unit 
vector in the direction of the magnetic field, and I is the unit 
tensor In equation (|3} an ideal Ohm's law is used, neglecting 
resistivity. In equation 0, P is the pressure tensor with dif- 
ferent perpendicular {p±) and parallel {p\\) components with 



respect to the background magnetic field, and 11 = 66(79 1| -p±) 
is the anisotropic stress tensor. (Note that 11 is not traceless in 
the convention used here.) P should in general be a sum over 
all species but in the limit where ion dynamics dominate and 
electrons just provide a neutralizing background, the pressure 
can be interpreted as the ion pressure. This is the case for hot 
accretion flows in which Tp ^ T^,. 

The exact pressures pn and p± can be rigorously de- 
termined by ta king moments of the drift kinetic equation 
jKulsrudlI983l) . 
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which is the asymptotic expansion of the Vlasov equation for 
the distribution function /^(x, /u, vy , f ) for species 's' with mass 
nis and charge in the limit p jL <^ 1, uj/fts ^ 1, where 
Ps and rts are the gyroradius and gyrofrequency, respectively. 
In equation (|5}, Ve = c(E x B) /B^ is the perpendicular drift 
velocity, /i = (v^ - Ve)^/2B is the magnetic moment (a con- 
served quantity in the absence of collisions), is the com- 
ponent of the gravitational force parallel to the direction of the 
magnetic field, and D/Dt = d/dt + (v^^b + \E) ■ V is the parti- 
cle Lagrangian derivative in configuration space. The fluid 
velocity V = Ve + 6m||, so the E x B drift is determined by the 
perpendicular component of equation Other drifts such 
as grad B, curvature, and gravity xB drifts are higher order 
in the drift kinetic MHD ordering and do not appear in this 
equation. In equation C(fs) is the collision operator to al- 
low for generalization to collisional regimes. Collisions can 
also be used to mimic rapid pitch angle scattering due to high 
frequency waves that break fi invariance. The parallel electric 

field is determined by = J^s^^f/ »Js)6 • V • Pj/ X]s("-5^s /'"O 
(Kulsrud 1983), which insures quasineutrality. 

Separate equations of state for the parallel and perpendic- 
ular pressures can also be obtained from the moments of the 
drift kinetic equation (Chew, Gold berger. & Lowlll956h . Ne- 
glecting the collision term these are: 



D 
pB— , 

Dt \pB, 

p3 D (p\\B^^ 



B2 Dt 



-V-q_L-^_LV-b, (6) 
-V-q||+2^^V-6, (7) 



where D/Dt = d/dt + Y ■ V is the fluid Lagrangian derivative 
and Qii^ = q\\_j_b are the heat fluxes (flux of and p±) par- 
allel to the magnetic field. The equation for the magnetic mo- 
ment density p{p) = p±/B can be written in a conservative 
form: 



If the heat fluxes are neglected (called the CGL or double adi- 
abatic limit), as the magnetic field strength (B) increases, p± 
increases (p± cx pB), and p|| decreases (p^^ cx p^ /B^). Inte- 
grating equation (jSJl over a finite periodic (even a shearing 
periodic) box shows that {p±/B) is conserved, where () de- 
notes a volume average. This implies that even when ^|| ^ 7^0, 
p± increases in a volume averaged sense as the magnetic en- 
ergy in the box increases. This means that that for a colli- 
sionless plasma, pressure anisotropy p± > {<)p\\ is created 
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as a natural consequence of processes that amplify (reduce) 
B. This pressure anisotropy is crucial for understanding mag- 
netic field amplification in collisionless dynamos. 

To solve the set of equations ( I1I4> . (I6I7> in a simple fluid 
based formalism, we require expressions for and q± in 
terms of lower order moments. No simple exact expressions 
for g'li and q± exist for nonlinear collisionless plasmas. Al- 
though simple, the double adiabatic or CGL approximation 
(where ^|| = = 0) does not capture key kinetic effects such 
as Landau damping. In the moderately collisional limit (p, < 
mean free path < system size), where the distribution func- 
tion is not very different from a local Ma xwellian, one can 
use the Braginskii equations (Braainskii 1965) to describe 
anisotropic transport (see Balbus 2000, 2004, for astrophysi- 
cal applications). However, in the hot RIAF regime, the mean 
free path is often much larger than the system size and the 
Braginskii equations are not formally applicable, though they 
are still useful as a qualitative indication of the importance of 
kinetic effects. The collisional limit of the kinetic MHD equa- 
tions can be shown to recover the dominant anisotropic heat 
flux and viscosity tensor of Braginskii (.SHD ). The local ap- 
proximation to kinetic MHD that we use here leads to equa- 
tions that are similar in form to Braginskii MHD, but with 
separate dynamical equations for parallel and perpendicular 
pressures. We also add models for enhanced pitch angle scat- 
tering by microinstabilities, which occur at very small scales 
and high frequencies beyond the range of validity of standard 
kinetic MHD. ' 

Hammett and collaborators have developed approximate 
fluid closures (c alled Lan dau fluid closure) for collisionless 
plasmas (H ammett & Perkins .1990: Hammett et al. 1 119921 
I^D) that capture kinetic effects such Landau damping. ISHDl 
give the resulting expressions for parallel heat fluxes (^||, ^j^) 
to be used in equations (|6j and Q. Landau closures are based 
on Pade approximations to the full kinetic plasma dispersion 
function that reproduce the correct asymptotic behavior in 
both the aj/^l|C|| 3> 1 and a'/A:||C|| <C 1 regimes (and provide a 
good approximation in between), where uj is the angular fre- 
quency, is the wavenumber parallel to the magnetic field, 

and C|| = y^p\\/p is the parallel thermal velocity of the parti- 
cles. In Fourier space, the linearized heat fluxes can be written 
as (equations (39) & (40) in SHD) 



^11 =■ 



POC\\Q 



-A)C||o 



'^11 {P\\/P) 
{p±/p) 



(9) 



+ 11 Ell [i-Ell^ 
^ " fio V P\\o . 



(10) 



where subscripts indicate equiUbrium quantities. Real space 
expressions are somewhat more cumbersome and are given by 
convolution integrals 



^11 



3/2 

- I nociio 



' This would also be needed when using Braginskii equations, because 
they are not necessarily well posed in situations where the anisotropic stress 
tensor can drive arbitrarily small scale instabilities! Schekochih in et al. 2005,). 
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' ^_, B(z + z')-B(z-z') 



f±0 



(12) 



where hq is the number density, = p\\/n and Tj_ = p±/n are 
the parallel and perpendicular temperatures, and z' is the spa- 
tial variable along the magnetic field line. SHO have shown 
that these fluid closures for the heat fluxes accurately repro- 
duce the kinetic linear Landau damping rate for all MHD 
modes (slow, Alfven, fast and entropy modes). The growth 
rate of the MRJ using the Landau closure model is also 
very similar to that obtained from full kinetic theory. As 
noted in the introduction, in addition to reproducing linear 
modes/instabilities. Landau fluid closures have also been used 
to model turbulence in fusion plasmas with reasonable suc- 
cess. 

These closure approximations were originally developed 
for turbulence problems in fusion energy devices with a strong 
guide magnetic field, where the parallel dynamics is essen- 
tially linear and FFTs could be easily used to quickly eval- 
uate the Fourier expressions above. In astrophysical prob- 
lems with larger amplitude fluctuations and tangled magnetic 
fields, evaluation of the heat fluxes become somewhat more 
complicated. One could evaluate the convolution expres- 
sions, equations (II 1> and ( I12l i (with some modest complex- 
ity involved in writing a subroutine to integrate along mag- 
netic field lines), leading to a code with a computational time 
Tcpii cx N^-N\\, where N^- is the number of spatial grid points 
and A^ii is the number of points kept in the integrals along field 
lines. (In some cases, it may be feasible to map the fluid quan- 
tities to and from a field-line following coordinate system so 
that FFTs can reduce this to Tcp,, oc N'} logA^n . ) While this is 
more expensive than simple MHD where Tcp,, cx N^., it could 
still represent a savings over a direct solution of the drift ki- 
netic equation, which would require Tcpu cx N^N^.^^ N^.^ , where 

Ny^^Ny^ is the number of grid points for velocity space. ^ 

As a first step for studying kinetic effects, in this paper we 
pick out a characteristic wavenumber ki that represents the 
scale of collisionless damping and use a local approximation 
for the heat fluxes in Fourier space, with a straightforward 
assumption about the nonlinear generalization: 



1\\=- 



8„„ V|i {pwIp) 

kL ' 

{p±/p) 



PC\ 



pc\\ 



-c\\P, 



1-^ 



V||B 



(13) 



(14) 



Note that this formulation of the heat flux is analogous to a 
Braginskii heat conduction along magnetic field lines. For lin- 
ear modes with | ^ ki, these approximations will of course 
agree with kinetic theory as well as the Pade approximations 

^ On the other hand, an effective hyperdiffusion operator in velocity space 
may reduce the velocity resolution requirements, and recent direct kinetic 
simulations of turbulence in fusion devices have found that often one does 
not need very high velocity resolution. This may make a direct solution of the 
drift kinetic equation tractable for some astrophysical kinetic MHD problems. 
Furthermore, a direct solution of the drift kinetic equation involves only local 
operations, and thus is somewhat easier to parallelize than the convolution 
integrals. 
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shown in One can think of ki as approximately con- 

trolUng the heat conduction rate, though this does not neces- 
sarily affect the resulting Landau damping rate of a mode in 
a monotonic way, since this sometimes exhibits impedance 
matching behavior. I.e., some modes are weakly damped in 
both the small and large (isothermal) heat conduction limits. 
We will vary to investigate the sensitivity of our results to 
this parameter. 

2.1. Linear Modes 

Since pressure anisotropy arises as a consequence of mag- 
netic field amplification in a collisionless plasma, it is of in- 
terest to repeat the linear analysis of the collisionless MRI 
done previously in .ODH,. but with a background pressure 
anisotropy {pp ^ p±o)- We consider the simple case of a ver- 
tical magnetic field. This analysis provides a useful guide to 
understanding some of our numerical results. 

We linearize equations for a differentially rotating 

disk (Vo = Rft(R)4)) with an anisotropic pressure about a uni- 
form subthermal vertical magnetic field (Bo = B^z). We as- 
sume that the background (unperturbed) plasma is described 
by a bi-Maxwellian distribution (pp ^ p±o)- We also as- 
sume that the perturbations are axisymmetric, of the form 
exp[-iujt + ik ■ x] with k = fegR + Lz. Writing p = po + Sp, 
B = Bo + (5B, p± = p±o + Sp±, p\\ = p||o + (5p||, working in 
cylindrical coordinates and making a \k\R ^ 1 assumption, 
the linearized versions of equations ([0-(|3} become: 

(15) 

SBR-ik^Sp^, (16) 

iP\\Q-p±o)\ 



LuSp= pok ■ S\, 
-iujpoSvR- polilSvcj, = 



ikR 



+ lk, { — - 
. 4-TT 



4tt 
iP\\o-p±o) 



-lujpodv^ + poOvr— = ik- I — ■ 



-iujp{)5v~ 

Lo5BR=—k^B^5vR 
uj6Bif,=-k,B^Sv^ 
LL>SB. = kRB^SvR, 



ikR{p\\o-pM))^-ikz5p\\, 



ik,B, dVL 
ijj dlnR 



Svr, 



5B^,{\1) 

(18) 
(19) 
(20) 
(21) 

where = 4^1^+ dil^ / dlnR is the epicyclic frequency. Equa- 
tions jl5>- J21> describe the linear modes of a collisionless 
disk with an initial pressure anisotropy about a vertical mag- 
netic field. This corresponds to the 9 = tt/2 case of ODH, 
but with an anisotropic initial pressure. Equations jl6> & 
Mil show that an initial anisotropic pressure modifies the 
Alfven wave characteristics, so we expect a background pres- 
sure anisotropy to have an important effect on the MRI. One 
way of interpreting equations M6\ & Ml\ is that p±^ > p\\ 
iP\\ > P-l) makes the magnetic fields more (less) stiff; as a 
result, this will shift the fastest growing MRI mode to larger 
(smaller) scales. 

The linearized equations for the parallel and perpendicu- 
lar pr essure response are given by equations (33) and (34) in 
ISHOt We present them here for the sake of completeness. 

-iu!5p\\ +p||o/k • 5\ + ik,q\\ + 2p\\(iik^5v, = 0, (22) 
— jw5/7j_ + 2/7_Lo!k • 5\+ik-qi_- p^Qik^Sv^ = 0, (23) 

where the heat fluxes can be expressed in terms of lower mo- 



ments using 



^c||o |^((5/?i -010(5/0)+ Y^C||o/?j_o X 

1 _;£M^ 

8" ik, _ 2 r 



(24) 
(25) 



where C|| q = \ /p\\ol po and SB = |(5B|. 

Figure lATI shows the MRI growth rate as a function of pres- 
sure anisotropy for two values of kR for (3 = 100. This figure 
shows that the fastest growing MHD mode {kR = 0) is stabi- 
lized for {pi_(i- P\\())/ P\\(i ^ 4//3; modes with kR^Q modes 
require larger anisotropy for stabilization. For /3 ^ 1, these 
results highlight that only a very small pressure anisotropy is 
required to stabilize the fastest growing MRI modes. Growth 
at large pressure anisotropics in Figure lAll for kR^Q mode 
is because of the mirror instability that is discussed below. 
The physical interpretation of the stabilization of the MRI in 
Figure lATl is that as the pressure anisotropy increases (p^o > 
P\\q), the field lines effectively become stiffer and modes of 
a given k can be stabilized (though longer wavelength modes 
will still be unstable). In a numerical simulation in which 
the pressure anisotropy is allowed (unphysically, as we see in 
§2.2) to grow without bound as the magnetic field grows, this 
effect is capable of stabilizing all of the MRI modes in the 
computational domain at very small amplitudes (see Fig. IA6I 
discussed in §4). 

2.2. Isotropization of the Pressure Tensor in Collisionless 
Plasmas 

Pressure anisotropy {p±^ ¥ P\\) is ^ source of free energy 
that can drive instabilities which act to isotropize the pressure, 
effectively providing an enhanced "co llision" rate in a colli- 
sionless plasma re.g.. lGarv et alJl997l) . In order to do so, the 
instabilities must break magnetic moment conservation, and 
thus must have frequencies comparable to the cyclotron fre- 
quency and/or parallel wavelengths comparable to the Larmor 
radius. Because of the large disparity in timescales between 
/i-breaking microinstabilities and the MRI (a;,„/„„/il ^ 10^), 
one can envision the microinstabilities as providing a "hard 
wall" limit on the pressure anisotropy: once the pressure 
anisotropy exceeds the threshold value where microinstabili- 
ties are driven and cause rapid pitch angle scattering, the pres- 
sure anisotropy nearly instantaneously reduces the anisotropy 
to its threshold value (from the point of view of the global 
disk dynamics). In this section we review the relevant insta- 
bilities that limit the pressure anisotropy in high (3 collision- 
less plasmas - these are the firehose, mirror, and ion cyclotron 
instabilities. We then discuss how we have implemented these 
estimated upper bounds on the pressure anisotropy in our nu- 
merical simulations. 

2.2.1. Maximum anisotropy for p\]^ > p± 

Plasmas with p\\ > p± can be unstable to the firehose in- 
stability, whose dispersion relation for p arallel propag ation is 
given by equation (2.12) of Kennel & Sagdee31l965: 



OJ 



u:n£.pj+nfkipf[i 



2;,2 2 



P\\ 



2 



= 0, 



(26) 



where = 87r/?||/B^, p,- is the ion Larmor radius, fli is the ion 
cyclotron frequency, and kn is the wavenumber parallel to the 
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local magnetic field direction. Solving for uj gives 



2 



± ik\\C\\ 



P±_ 
P\\ 



2 



Ir ' 



1/2 



(27) 



For long wavelengths, the firehose instability requires /?|| > 
P±+B^/4tt and is essentially an Alfven wave destabilized 
by the pressure anisotropy. The maximum growth rate oc- 
curs when k^^pj = 2{1- p±/ pn -2//3||) and is given by 17,(1 ~ 
px/pil - 2//9||). We use an upper limit on pn > p± corre- 
sponding to l-/7x//'||-2//3|| < 1/2, which is an approximate 
condition for the growth of modes that will violate fi conser- 
vation and produce rapid pitch angle scattering. 



2.2.2. Maximum anisotropy for p± > p 



For p± > there are two instabilities that act to isotropize 
the pressure, the mirror instabi lity and the ion cyclotron in- 
stability ('e.g.. lGarv et alJ[T997h . A plasma is unstable to the 
mirror instability when p±/p\\ - 1 > 1 /(3±, although as dis- 
cussed below only for somewhat larger anisotropies is mag- 
netic moment conservation violated. Formally, a plasma with 
any nonzero press ure an isotropy can be unstable to the ion cy- 
clotron instability (lStixlfr99 2). However, there is an effective 



threshold given by the requirement that the unstable modes 
grow on a timescale comparabl e to the disk r otati on period. 

Equations (43') & (44') of iHasegawal /T969) give the 
wavenumber for the fastest growing mirror mode 



k±pi- 



(D-l) 



(D-l) 



(28) 
(29) 



where D = (3±ip±/p\\ - 1), f3± = Sttp±/B^. To estimate the 
pressure anisotropy at which fi conservation is broken and 
thus pitch angle scattering is efficient, we calculate D for 
which ^ k±pi ^ 1. This implies D w 7 or that p, conser- 
vation fails (and pitch angle scattering occurs) if the pressure 
anisotropy satisfies 



p\\ Pi- 



rn 



The ion cyclotron instability can be also be excited when 
pi. > P\\- Gary and collaborators have analyzed the ion cy- 
clotron instability in detail through l inear analysis and numer- 
ical simulations. iGarv et al.1 \\991\i and iGarv & Led J 1994 
calculate the pressure anisotropy required for a given growth 
rate 7 relative to the ion cyclotron frequency i7, 



P\\ 



1 > 



PI 



(31) 



where S' = 0.35 and p = 0.42 are fit ting parameters quoted in 
equation (2) of lGarv&Le^ jl994 for 7/r2,- =10"'*. More- 
over, for 7 <C ri, the threshold anisotropy depends only very 
weakly on the growth rate 7. As a result, equation ( 13 U pro- 
vides a reasonable estimate of the pressure anisotropy re- 
quired for pitch angle scattering by the ion cyclotron insta- 
bility to be important on a timescale comparable to the disk 
rotation period. 



2.3. Pressure Anisotropy Limits 

Motivated by the above considerations, we require that the 
pressure anisotropy satisfy the following inequalities in our 
simulations (at each grid point and for all time steps): 



P\\ 2' 
P 



P\\ ^ 



P\\ 



1/2 



(32) 
(33) 

(34) 



where S and ^ are constants described below. It is important 
to note that the fluid-based kinetic theory utilized in this paper 
can correctly reproduce the existence and growth rates of the 
firehose and mirror instabilities (though not the ion cyclotron 
instability).^ However, it can only do so for long wavelength 
perturbations that conserve p. The relevant modes for pitch 
angle scattering occur at the Larmor radius scale, which is 
very small in typical accretion flows and is unresolved in 
our simulations. For this reason we must impose limits on 
the pressure anisotropy and cannot simultaneously simulate 
the MRI and the relevant instabilities that limit the pressure 
anisotropy. The algorithm to i mpos e the pressure anisotropy 
limits is explained in ApDendix lA.3l 

In equation j33> . the parameter ^ determines the threshold 
anisotropy above which the mirror instability leads to pitch 
angle scattering. A value of ^ = 3.5 was estimated in § 12.2.21 
We take this as our fiducial value but for comparison also 
describe calculations with = 0.5, which corresponds to the 
marginal state for the mirror instability. We compare both 
models because the saturation of the mirror instability is not 
well understood, particularly under the conditions appropri- 
ate to a turbulent accretion disk. Equation j34> i s based on 
the pitch angle scattering model used by Birn & Hessd J2001I) 
for simulations of magnetic reconnection in collisionless plas- 
mas; following them we choose 5 = 0.3. Equation J34t with 
5 = 0.3 gives results which are nearly identical (for the typical 
range of (3 studied here) to the pressure anisotropy threshold 
for the ion cyclotron instability discussed in ^2.2.2l (eq. l31l ). 

In our simulations we find that for typical calculations if 
^ = 0.5 then equation ( I33t (the "mirror instability") domi- 
nates the isotropization of the pressure tensor while if ^ = 3.5 
then equation j34t (the "ion cyclotron instability") dominates. 
We also find that our results are insensitive to the form of 
the p\\ > p± threshold (eq. 1321 ): e.g., simulations with 
l-p±/p\\ < 2//3[| (the marginal state of the firehose mode) 
instead of equation (j^D give nearly identical results. Future 
fully kinetic simulations of the mirror, firehose, and ion cy- 
clotron instabilities will be useful for calibrating the pitch an- 
gle scattering models used here. 

3. NUMERICAL METHODS 

In this section we discuss the shearing box equations that 
we solve numerically and the modifications made to ZEUS to 
include kinetic effects. 

3.1. Shearing Box 



The double adiabatic limit {q± = =0) predicts an incorrect threshold 
and incorrect growth rates for the mirror instability (e.g., SHD). Thus it is 
important to use the heat flux models described in §2 to capture the physics 
of the mirror instability. 
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The shearing box is based on a local expansion of the tidal 
forces in a reference frame corotating with the disk (see HGB 
for details). A fiducial radius TJq in the disk is picked out and 
the analysis is restricted to a local Cartesian patch such that 
Lx,Ly,L~ <C Ro (where x = r-Ro, y = (j) and z = z). In this 
paper only the radial component of gravity is considered and 
buoyancy effects are ignored. We also assume a Keplerian 
rotation profile. With these approximations, the equations of 
Landau MHD in the shearing box are: 

dp 

'^. + V-(pV) = 0, 



dt 

dY 1 / B-VB 

— +V- VV = — V I p± + — + 

--\/ ■Il-2nxY + 3n^xx, 
P 

9B 

— = Vx(VxB), 



(35) 

(36) 
(37) 



^ + V • (pii V) + V • qii +2p||b • VV • h-lq^y ■ b 

dt " " " 



dp± 
dt 



= -:^^eff(p\\-p±), 
+ V • (pi_Y) + V ■ qj_ + p^V -Y-pAh • VV • b 

+ q±^ ■b = -'^i^eff(p±-p\\), 



^_l=-Pk_lV[| ( )+K„,B-VB, 



(38) 

(39) 
(40) 

(41) 



where q|| = ^||b and = q±b are the heat fluxes parallel to 
the magnetic field, Veff is the effective pitch-angle scatter- 
ing rate (including microinstabilities, see §2.3 and Appendix 
IA.3> . K|| and k± are the coefficients of heat conduction, and 
K,„ is the coefficient in q± due to parallel gradients in the 
strength of magnetic field (SHD). The k,,, component of 
q± that arises because of parallel magnetic field gradients is 
important for correctly recovering the saturated state for the 
mirror instability in the fluid limit, where (in steady state) 
^11 X w implies that is constant along the field line, and 
T± and magnetic pressure are anticorrelated. 

Given our closure models, the coefficients for the heat flux 
are given by 



=11 



1 



87:fkL + {3n-8)iy,ff 



K± = 



1 \ P^ 



(42) 
(43) 
(44) 



where ki is the parameter that corresponds to a typical 
wavenumber characterizing Landau damping (see §2). We 
consider several values of ki to study the effect of Lan- 
dau damping on different scales. In particular, we consider 
kL = 0.5 / &, 0.25 /5z, 0. 125/ 5z which correspond to correctly 
capturing Landau damping on scales of 12&, 24&, 48&, re- 
spectively, where 5z = L./N^, L^ = l for all our runs, and A^, is 



the number of grid points used in the z-direction (taken be 27 
and 54 for low and high resolution calculations, respectively). 
Thus, ki = 0.25/ 6z corresponds to correctly capturing Landau 
damping for modes with wavelengths comparable to the size 
of the box in the low resolution runs. 

The term i^gjf in equations J42> & j43> is an effective colli- 
sion frequency which is equal to the real collision frequency 
ly as long as p conservation is satisfied. However, when the 
pressure anisotropy is large enough to drive microinstabilities 
that break /i invariance and enhance pitch angle scattering, 
then there is an increase in the effective collision frequency 
that decreases the associated con ductiv i ties. T he ex pressi ons 
for i^fff a re gi ven in equations (IA12> . ( IA15> , and jA18> of 
Appendix lA.3l 

Shearing periodic boundary conditions appropriate to the 
shearing box are described in HGB. Excluding V,, all vari- 
ables at the inner x- boundary are mapped to sheared ghost 
zones at the outer boundary; a similar procedure applies for 
the inner ghost zones. V,. has a jump of (3/2)i7Li across the 
box while applying the x- shearing boundary conditions, to 
account for the background shear in Vy. 

3.2. Numerical Methods 

We have used a version of the ZEUS code modified to in- 
clude kinetic effects (see Stone & Norman 1992a b). ZEUS 
is a time explicit, operator split, finite difference algorithm 
on a staggered mesh, i.e., scalars and the diagonal compo- 
nents of second rank tensors are zone centered, while vectors 
are located at zone faces, and pseudovectors and offdiagonal 
components of second rank tensors are located at the edges. 
The location of different variables on the grid is described in 
more detail in Appendix lA.il Appendix IA.2I describes how 
we choose the time step 6t to satisfy the Courant condition 
(which is modified by pressure anisotropy and heat conduc- 
tion). We also require that the choice of 6t maintain positivity 
of and p±. 

Implementation of the shearing box boundary conditions is 
described in HGB. One can either apply boundary conditions 
on the components of B or the EMF. We apply shearing peri- 
odic boundary conditions on the EMF to preserve the net ver- 
tical flux in the box, although applying boundary conditions 
directly on B also gives satisfactory results. 

Equations ( I38> and ( I39l l are split into a transport and source 
step, analogous to the energy equation in the original MHD 
formalism. The transport step is advanced conservatively, and 
source step uses central differences in space. It should be 
noted that in equation ( I39> the V • q^ term is not purely diffu- 
sive, and it is necessary to carefully treat the magnetic gradi- 
ent part of q± in the transport step for robustness of the code 
(Appendix lA.4> . 

We have carried out a series of tests of our newly added 
subroutines for evolving anisotropic pressure and parallel 
heat conduction. We tested the anisotropic conduction rou- 
tine by initializing a "hot" patch in circular magnetic field 
lines and assessing the extent to which heat remains con- 
fined along the field. This is the same test described in de- 
tail in Parrish & Stone (2005) and we find good agreement 
with their results. For the diffusion of a narrow temperature 
pulse in ID, the anisotropic conduction routine also gave re- 
sults nearly identical to that predicted analytically by the ID 
diffusion equation.** Additional tests of the code included lin- 

F or details of the test and error analysis see: 
IhttpV/wS.pppl.gov/'psharma/cartesian/ lDheatdit'fusion/ 
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ear waves and instabilities in non-rotating anisotropic plas- 
mas, including the Alfven wave and the firehose and mirror 
instabilities. For mirror simulations we observe the formation 
of stationary anticorrelated density and magnetic structures as 
seen i n the hybrid simulations of McKean, Garv, & Winske 
^199%. For firehose we see the instability with magnetic per- 
turbations developing at small scales but during satur ation the 
perturb ations are at larger scales as seen in lOuest & Shapirol 

Finally, the numerical growth rates of the kinetic MRI 
were compared to the analytic results for different pressure 
anisotropics, (kx,kz), collision frequencies, and angles be- 
tween the magnetic field and z; we find good agreement with 
the results of ODH and SHO. When = k^^, the growth rate 
of the fastest growing mode is within 3% of the theoreti- 
cal prediction. The case Z?^ = B. shows ~ 2 faster growth as 
compared to = as predicted by linear theory. 

3.3. Shearing Box and Kinetic MHD 

Certain analytic constraints on the properties and energet- 
ics of shearing box simulations have been described in iHGBl 
These constraints serve as a useful check on the numerical 
simulations. Here we mention the modifications to these con- 
straints in KMHD. Conservation of total energy in the shear- 
ing box gives 



dt 2 



dA 



pV.SV,- 



47r(/?||-pi)\ B,By 



B2 



4tt 



(45) 

where dVy = Vv + (3 /2)ilx, and F is the total energy given by. 



r = 



d^x 



P 



2 ^ 



(46) 



where (/> = -3/217^ is the tidal effective potential about /Jo- 
Equation ( 145 I I states that the change in the total energy of 
the shearing box is due to work done on the box by the 
boundaries. Notice that there is an anisotropic pressure con- 
tribution to the work d one on the box. Equation (29) in 
iBalbus & Hawlevljl998h for conservation of angular momen- 
tum in cyUndrical geometry is also modified because of the 
anisotropic pressure and is given by 



^(p/?y^)+v. 



4ir{p\\-pi_) 
B2 



= 0, 



(47) 



where Bp = BrR + Z?;Z is the poloidal field. We can calculate 
the level of angular momentum transport in our simulations 
by measuring the stress tensor given by 

w„ = py.5y„-Mi + tei)5A (48) 

Ait B^ 

Note that the stress tensor has an additional contribution due 
to pressure anisotropy. One can define a dimensionless stress 
via Shakura and Sunyaev's a parameter by 

«=^^ 



■ aR + am + 



(49) 



where aR, olm, ctA are the Reynolds, Maxwell and anisotropic 
stress parameters, respectively. As in iHGBl we normalize the 
stress using the initial pressure to define an a parameter 



3.4. Shearing Box Parameters and Initial Conditions 

The parameters for our baseline case have been chosen to 
match the fiducial run Z4 of HGB. The simulation box has 
a radial size L^ = 1, azimuthal size Ly = 2tt, and vertical size 

L, = 1 . The sound speed y = \/p/p = LSI, so that the vertical 
size is about a disk scale height (though it is an unstratified 
box). The pressure is assumed to be isotropic initially, with 
= pqV} = 10"^ and po = 1- All of our simulations start with 
a vertical field with j3 = Stt/jo/Bq = 400. The fastest growing 
MRI mode for this choice of parameters is reasonably well 
resolved. We consider two different numerical resolutions: 
27 X 59 X 27, and 54 x 1 18 x 54. Perturbations are introduced 
as initially uncorrelated velocity fluctuations. These fluctua- 
tions are randomly and uniformly distributed throughout the 
box. They have a mean amplitude oi\5V\ = 10"^,. 

4. RESULTS 

The i mpor tant parameters for our simulations are listed 
in Table lAll Each simulation is labeled by Z (for the ini- 
tial field), and I and h represent low (27 x 59 x 27) and 
high (54 X 118 X 54) resolution runs, respectively. We also 
include low and high resolution MHD runs for comparison 
with the kinetic calculations (labeled by ZM). Our models 
for heat conduction and pressure isotropization have several 
parameters: ki, the typical wavenumber for Landau damp- 
ing used in the heat flux (eqs. 1131 & 1141 ). and ^, the pa- 
rameter that forces the pressure anisotropy to be limited by 
P^l fll ~ 1 < "^^1 Pi- (representing pitch angle scattering due 
to small scale mirror modes; eq. 1331 ). All of our calculations 
except Z/8, Z/1, and Zh\ also use the ion cyclotron scattering 
"hard wall" from equation ( I34t . In addition to these model 
parameters. Table lAll also lists the results of the simulations, 
including the volume and time averaged magnetic and kinetic 
energies, and Maxwell, Reynolds, and anisotropic stresses. 
As Table lAll indicates. the results of our simulations depend 
quantitatively - though generally not qualitatively - on the 
microphysics associated with heat conduction and pressure 
isotropization. Throughout this section we use single brackets 
(/) to denote a volume average of quantity /; we use double 
brackets ((/)) to denote a volume and time average in the sat- 
urated turbulent state, from orbit 5 onwards. 

4.1. Fiducial Run 

We have selected run Z/4 as our fiducial model to describe 
in detail. This model includes isotropization by ion cyclotron 
instabilities and mirror modes, with the former dominating 
(for ^ = 3.5; see ^2.2.21 except at early times. The conduc- 
tivity is determined by k^ = 0.5/& which implies that modes 
with wavelengths ~ 12& ^ L,/2 are damped at a rate consis- 
tent with linear th eory. 

Fi gure s lA2l A4I show the time evolution of various physical 
quantities for run ZZ4. The early linear development of the 
instability is similar to that in MHD, with the field growing 
exponentially in time. The key new feature is the simultane- 
ous exponential growth of pressure anisotropy (p^ > p \\) as 
a result of p conservation (up to 2 orbits in Fig. IA4> . As 
described in ^2.11 this pressure anisotropy tends to stabilize 
the MRI modes and shut off the growth of the magnetic field. 
Indeed, in simulations that do not include any isotropization 
of the pressure tensor, we find that all MRI modes in the box 
are stabilized by the pressure anisotropy and the simulation 
saturates with the box filled with small amplitude anisotropic 
Alfven waves (see Fig. lA6> . This highlights the fact that, un- 
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like in MHD, the MRI is not an exact nonlinear solution in 
kinetic theory. However, the pressure anisotropy required to 
stabilize all MRI modes exceeds the pressure anisotropy at 
which pitch angle scattering due to mirror and ion cyclotron 
instabilities become important. This takes place at about orbit 
2 in run Z/4 (see the small 'dip' in the growth of B in Fig. IA2> . 
at which point the pressure anisotropy is significantly reduced 
and the magnetic field is able to grow to nonlinear amplitudes. 

The nonlinear saturation at orbit ~ 5 appears qualitatively 
similar to that in MHD, and may o ccur via analogues of th e 
parasitic instabilities described by iGoodman & Xul Jl994h . 
The channel solution is, however, much m ore extreme in 
KMHD than MHD (the maximum in Fig. IA2l is approxi- 
mately an order of magnitude larger than in analogous MHD 
runs). After saturation, the magnetic and kinetic energies in 
the saturated state are comparable in KMHD and MHD (see 
Table lAll . This is essentially because the pitch angle scatter- 
ing induced by the kinetic microinstabilities acts to isotropize 
the pressure, enforcing a degree of MHD-like dynamics on 
the collisionless plasma^ 

Figure |A3] and Table lAll show the various contributions to 
the total stress. As in MHD, the Reynolds stress is signifi- 
cantly smaller than the Maxwell stress. In kinetic theory, how- 
ever, there is an additional component to the stress due to the 
anisotropic pressure (eq. 1471 ). In the saturated state, we find 
that the Maxwell stress is similar in KMHD and MHD, but 
that the anisotropic stress itself is comparable to the Maxwell 
stress. Expressed in terms of an a normalized to the initial 
pressure, our fiducial run Z/4 has = 0.23, = 0.097, and 
ttA = 0.2, indicating that stress due to pressure anisotropy is 
dynamically important. 

Nearly all physical quantities in Figures lA2IA4l reach an ap- 
proximate statistical steady state. The exceptions are that p\\ 
and p±^ increase steadily in time because the momentum flux 
on the boundaries does work on the system (Eq. 45), which is 
eventually converted to heat in the plasma by artificial viscos- 
ity and there is no cooling (the same is true in HGB's MHD 
simulations). Because of the steadily increasing internal en- 
ergy and approximately fixed (although with large fluc- 
tuations), the plasma (3 shows a small secular increase from 
orbits 5-20 (a factor of w 3 increase, though with very large 
fluctuations due to the large fluctuations in magnetic energy). 
Figure lA4l shows the pressure anisotropy thresholds due to the 
ion cyclotron and mirror instabilities, in addition to the vol- 
ume averaged pressure anisotropy in run Z/4. From equation 
( I34> . the ion cyclotron threshold is expected to scale as 
which is reasonably consistent with the trend in Figure lA4l 
The actual pressure anisotropy in the simulation shows a small 
increase in time as well, although less than that of the ion cy- 
clotron threshold. These secular changes in /3 and Ap are a 
consequence of the increasing internal energy in the shear- 
ing box, and are probably not realistic. In a global disk, we 
expect that - except perhaps near the inner and outer bound- 
aries - (3 will not undergo significant secular changes in time. 
In a small region of a real disk in statistical equilibrium, the 
heating would be balanced by radiation or by cooler plasma 
entering at large R and hotter plasma leavi ng at small R. 

It is interesting to note that in Figure IA4I the pressure 
anisotropy (A-nAp/B^) is closely tied to the ion cyclotron 
threshold at times when B^ is rising (which corresponds to the 
channel solution reemerging). Increasing B leads to a pressure 
anisotropy with p± > p^^ by p conservation. At the same time, 
the ion cyclotron threshold (~ y/J5) decreases and thus the 



threshold is encountered which limits the anisotropy. When 
B is decreasing, however, we do not find the same tight re- 
lationship between the pressure anisotropy and the imposed 
threshold. Figure IA4I clearly indicates that in our fiducial 
simulation pitch angle scattering is dominated by the ion cy- 
clotron threshold. For comparison, Figure lX51 shows the pres- 
sure anisotropy and thresholds for run Z/8 which is identical 
to the fiducial run except that the ion cyclotron threshold is not 
used and the only scattering is due to the mirror threshold. In 
this case, the saturated pressure anisotropy is somewhat larger 
than in the fiducial run, but the pressure anisotropy is not tied 
to the mirror threshold. 

Table IA2l gives the mean, standard deviation, and standard 
error in the mean, for various quantities in the saturated por- 
tion of the fiducial simulation. The standard errors are esti- 
mated by taking into account the finite correlation time for 
the physical quantities in the simulation, as described in Ap- 
pendix IA.5I In many cases, the deviations are significantly 
larger than the mean. As in MHD, we find that the magnetic 
energy is dominated by the y-component, which is about a 
factor of 3 larger than the jc-component; the vertical compo- 
nent is smaller yet. The radial and azimuthal kinetic energy 
fluctuations are comparable, while the vertical component is 
smaller We also find that, as in MHD, the perturbed kinetic 
and magnetic energies are not in exact equipartition: the mag- 
netic energy is consistently larger Table lA2l also shows the 
mean and deviations for {p±/B) and (p^^B^/p^). Because 
of pitch angle scattering /i = {p±/B) is no longer conserved. 
{p\\B^ / p^) varies both because of heat conduction and pitch 
angle scattering. 

The pressure anisotropy in our fiducial run saturates at 
A--K{p±-p\\)/B^ K, 1.5. By contrast, the threshold for the mir- 
ror instability is 4Tr{p±- p0/B^ = 0.5. This implies that the 
model is unstable to generating mirror modes. However, the 
mirror modes that can be excited at this level of anisotropy do 
not violate p conservation and thus do not contribute to pitch 
angle scattering ( ^2.2.2> . They can in principle isotropize the 
plasma in a volume averaged sense by spatially redistribut- 
ing pl asma into magnetic wells (e.g.. lKivelson & South woodi 
This saturation mechanism can be calculated using our 
kinetic-MHD code and was in fact one of our test problems 
(for a uniform plasma). It does not appear to be fully efficient 
in the saturated state of our turbulent disk simulations, even at 
the highest resolutions we have run. 

In the next few sections we compare the fiducial simulation 
described above with variations in the pitch angle scattering 
model and the parallel conductivity. A compariso n of the total 
stress in all of our simulations is shown in Figure lXTl 

4.2. The Double Adiabatic Limit 

Simulations Z/1 and Z/il are simulations in the double adi- 
abatic limit (no heat conduction), with no limit on the pres- 
sure anisotropy imposed. In this limit both /i = {p±/B) and 
{p^B^ / p^) are conserved. Figure lXSl shows volume averages 
of various quantities as a function of time for the run Z/1. 
These calculations are very different from the rest of our re- 
sults and show saturation at very low amplitudes {5B^ /B^ w 
0.04). In the saturated state, the box is filled with shear mod- 
ified anisotropic Alfven waves and all physical quantities are 
oscillating in time. The total stress is also oscillatory with 
a vanishing mean, resulting in negligible transport. In these 
calculations, the pressure anisotropy grows to such a large 
value that it shuts off the growth of all of the resolved MRI 
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modes in the box. Table lAll shows that ((47r(/?|| - p±)/B^)) 
saturates at -11.96 and -10.2 for the low and high reso- 
lution runs, respectively (although the normalized pressure 
anisotropy {{(p\\- p±)/p\\)) w -0.07 is quite small). This 
is much larger than the anisotropy thresholds for pitch angle 
scattering described in §2. As a result, we do not expect these 
cases to be representative of the actual physics of collisionless 
disks. These cases are of interest, however, in supporting the 
predictions of the linear theory with anisotropic initial condi- 
tions considered in ii l2.1l and in providing a simple test for the 
simulations. They also highlight the central role of pressure 
isotropization in collisionless dynamos. 

4.3. Varying Conductivity 

We have carried out a series of simulations with different 
conductivities defined by the parameter k^. Simulations Z/2 
and Zh2 are in the CGL limit with vanishing parallel heat con- 
duction, but with the same limits on pressure anisotropy as 
the fiducial model. Simulations Z6 use kiSz = 0.25 while run 
Zn uses ki = 0.125/&. Both of these are smaller than the 
value of kiSz = 0.5 i n the fiducial run, which implies a larger 
conductivity. Figure IATI shows that the total stress varies by 
about a factor of 2 depending on the conductivity and res- 
olution. Simulations with larger conductivity tend to have 
smaller saturation amplitudes and stresses. This could be be- 
cause larger conductivity implies more rapid Landau damping 
of slow and fast magnetosonic waves. In all cases, however, 
the anisotropic stress is comparable to the Maxwell stress as in 
the fiducial run. Until a more accurate evaluation is available 
of the heat fluxes for modes of all wavelengths in the simu- 
lation simultaneously (either by a more complete evaluation 
of the nonlocal heat fluxes, eqs. II 11121 . or even by a fully 
kinetic MHD code that directly solves eq. ||5l, it is difficult to 
ascertain which value of the conductivity best reflects the true 
physics of collisionless disks. 

4.4. Different Pitch Angle Scattering Models 

In this section we consider variations in our model for pitch 
angle scattering by high frequency waves. All of these calcu- 
lations utilize ki = 0.5 /6z- We note again that the appropri- 
ate pitch angle scattering model remains somewhat uncertain, 
primarily because of uncertainties in the nonlinear saturation 
of long-wavelength /i-conserving mirror modes. The calcula- 
tions reported here cover what we believe is a plausible range 
of models. 

Models Z15 and Zh5 place a more stringent limit on the 
allowed pressure anisotropy, taking ^ = 0.5 in equation (I33> . 
This corresponds to the threshold of the mirror instability. Not 
surprisingly, this simulation is the most "MHD-like" of our 
calculations, with magnetic and kinetic energies and Maxwell 
stresses that are quite similar to those in MHD. Even with 
this stringent limit, however, the anisotropic stress is « 1/3 
of the Maxwell stress. It is also interesting to note that al- 
though the dimensionless pressure anisotropy is quite small 
{{4tt(p\^- p±)/B^)) w -0.02, the dimensionless anisotropic 
stress ((47r(p|| -p±)/B^ x B^By/po)) « -0.07 is significantly 
larger (and larger than Reynolds stress) because of correla- 
tions between the pressure anisotropy and field strength. 

As a test of how large a collisionality is needed for the 
results of our kinetic simulations to rigorously approach the 
MHD limit, we have carried out a series of simulations in- 
cluding an explicit collisionality v and varying its magnitude 
relative to the disk frequency il. Our results are summarized 



in Table lA3l and Figure lASl In these simulations we start with 
initial conditions determined by the saturated turbulent state 
of our fiducial run Z/4, but with an explicit collision freq uency 
(in addition to the scattering models described in ^2.3> . Fig- 
ure |A8l shows that for v/Vt < 20, the results are very simi- 
lar to the collisionless limit. For larger collision frequencies 
the anisotropic stress is reduced and the simulations quantita- 
tively approach the MHD limit. These results are similar to 
those obtained by SHO, who found that in linear calculations 
the MHD limit for modes with A: f2/vyi is approached when 

To consider the opposite limit of lower collisionality, run 
Z/8 places a less stringent limit on the allowed pressure 
anisotropy, taking ^ = 3.5 in equation ( I33> and ignoring the 
limit set by the ion cyclotron instability in equation ( I34> . The 
results of this calculation are not physical but are useful for 
further clarifying the relative importance of the Maxwell and 
anisotropic stresses as a function of the pitch angle scatter- 
ing rate. In Z/8, the saturated magnetic energy and Maxwell 
stress are lower than in all of our other calculations (excluding 
the double adiabatic models described in ^4.2> . Interestingly, 
however, the tot al str ess is comparable to that in the other cal- 
culations (Fig. IA7> because the anisotropic stress is w 2.4 
times large r tha n the Maxwell stress (Table 1). As discussed 
briefly in § 14.11 the pressure anisotropy in this simulation is 
not simply set by t he a pplied mirror pitch angle scattering 
threshold (see Fig. IA5t . It is possible that resolved mirror 
modes contribute to decreasing the volume averaged pressure 
anisotropy (but see below). 

Finally, in models Z3 we include parallel heat conduction 
but do not limit the pressure anisotropy. In these calcula- 
tions, we expect to be able to resolve the long-wavelength /i- 
conserving mirror modes that reduce the pressure anisotropy 
by fo rming magnetic wells (as in Kivelson & Southwoo^ 
119961) .^ In our test problems with uniform anisotropic plas- 
mas, this is precisely what we find. In the shearing box calcu- 
lations, however, even at the highest resolutions, we find that 
the pressure anisotropy becomes so large that equations ( I33> 
and ( I34> are violated and pitch angle scattering due to high 
frequency microinstabilities would become important. The 
resolved mirror modes are thus not able to isotropize the pres- 
sure sufficiently fast at all places in the box.^ However, it is 
hard to draw any firm conclusions from these simulations be- 
cause they stop at around 4 orbits (for both resolutions ZZ3 
and ZhJ) during the initial nonlinear transient stage. At this 
time the pressure becomes highly anisotropic and becomes 
very small at some grid points, and the time step limit causes 
5t~^Q. 

5. SUMMARY & DISCUSSION 

In this paper we have carried out local shearing box sim- 
ulations of the magnetorotational instability in a collisionless 
plasma. We are motivated by the application to hot radiatively 
inefficient flows which are believed to be present in many low- 
luminosity accreting systems. Our method for simulating the 
dynamics of a collisionless plasma is fluid-based and relies 
on evolving a pressure tensor with closure models for the heat 
flux along magnetic field lines (§2). These heat flux mod- 

^ At the resolution of Z/3, the fastest growing mirror mode in the compu- 
tational domain has a linear growth comparable to that of the MRI. 

* In higher resolution simulations, one can resolve smaller-scale and faster 
growing miiTor modes, and thus the effects of isotropization by resolved mir- 
ror modes could be come increasingly important. We see no such indications, 
however, for the range of resolutions we have been able to simulate. 
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els can also be thought of as approximating the collisionless 
(Landau) damping of linear modes in the simulation. 

By adiabatic invariance, a slow increase (decrease) in 
the magnetic field strength tends to give rise to a pressure 
anisotropy with p± > (py > p±), where the directions 
are defined by the local magnetic field. Such a pressure 
anisotropy can, however, give rise to small scale kinetic in- 
stabilities (firehose, mirror, and ion cyclotron) which act to 
isotropize the pressure tensor, effectively providing an en- 
hanced rate of pitch angle scattering ("collisions"). We 
have included the effects of this isotropization via a subgrid 
model which restricts the allowed magnitude of the pressure 
anisotropy ( ^2.3> . 

We find that the nonlinear evolution of the MRI in a col- 
lisionless plasma is qualitatively similar to that in MHD, 
with comparable saturation magnetic field strengths and mag- 
netic stresses. The primary new effect in kinetic theory 
is the existence of angular momentum transport due to the 
anisotropic pressure stress (eq. 1471 ). For the allowed pres- 
sure anisotropics estimated in ^2.31 the anisotropic stress is 
dynamically important and is as large as the Maxwell stress 
(Table 1). The precise rate of transport in the present simu- 
lations is difficult to quantify accurately and depends - at the 
factor of ^ 2 level - on some of the uncertain microphysics 
in our kinetic analysis (e.g., the rate of heat conduction along 
magnetic field lines and the exact threshold for pitch angle 
scattering by small-scale instabilities; see Fig. |^}. For 
better results, it would be interesting to extend these calcu- 
lations with a more accurate evaluation of the actual nonlo- 
cal heat fluxes, equations Jll> - J12> . or even to directly solve 
equation (|5j for the particle distribution function. Further ki- 
netic studies in the local shearing box, including studies of the 
small-scale instabilities that limit pressure anisotropy, would 
be helpful in developing appropriate fluid closures for global 
simulations. 

It is interesting to note that two-temperature RlAFs can 
only be maintai ned below a critical luminosity ^ o^Ledd 
SRses et. alJl982 ). Thus enhanced transport in kinetic theory 
due to the anisotropic pressure stress would extend upward in 
luminosity the range of systems to which RIAFs could be ap- 
plicable. This is important for understanding, e.g., state tran- 
sition s in X-ray binaries (e.g.. lEsin. McClintock. & NaravanI 

In addition to angular momentum transport by anisotropic 
pressure stresses. Landau damping of long-wavelength modes 
can be dynamically important in collisionless accretion flows. 
Because the ZEUS code we employ is non-conservative, we 
cannot carry out a rigorous calculation of heating by different 
mechanisms such as Landau damping and reconn ection. Fol- 
lowing the total energy-conserving scheme of Turn er et'aP 
(|2003), however, we estimate that the energy dissipated by 
Landau damping is comparable to or larger than that due re- 
connection (which is the major source of heating in MHD 
simulations). One caveat to this analysis is that in local sim- 
ulations, the pressure increases in time due to heating, while 
^ constant. Thus /3 increases in time and the turbulence 
becomes more and more incompressible. This will artificially 



decrease the importance of compressible channels of heating. 
Clearly it is of significant interest to better understand heat- 
ing and energy dissipation in RIAFs, particularly for the elec- 
trons. We will carry out a more a systematic analysis of the 
energetics of collisionless disks in future global simulations. 

In all of our calculations, we have assumed that the dom- 
inant source of pitch angle scattering is high frequency mi- 
croinstabilities generated during the growth and nonlinear 
evolution of the MRI. We cannot, however, rule out that there 
are other sources of high frequency waves that pitch angle 
scatter and effectively decrease the mean free path of parti- 
cles relative to that calculate d he re (e.g., sh ocks and recon- 
nection). As shown in Table IX^ and Figure \KE\ this would 
decrease the magnitude of the anisotropic stress; we find that 
for > 30 57, the results of our kinetic simulations quantita- 
tively approach the MHD limit. In this context it is impor- 
tant to note that the incompressible part of the MHD cascade 
launched by the MRI is expected to be highly anisotropic with 
> k\\ (Goldreich& Sridhar 1995). As a result, there is 
very little power in high frequency waves that could break /i 
conservation. It is also interesting to note that satellites have 
observed that the pressure anisotropy in the solar wind near 1 
AU is approximately marginally stable to the firehose instabil- 
ity ( KasDcr et al. 2002), consistent with our assumption that 
microinstabilities dominate the isotropization of the plasma. 

In this paper we have focused on kinetic modifications to 
angular momentum transport via anisotropic pressure stresses 
and parallel heat conduction. In addition, kinetic effects sub- 
stantially modify the stability of thermally stratified low col- 
lisionality plasmas such as those expected in RIAFs. Balbi^ 
(,2000) showed that in the presence of anisotropic heat con- 
duction, thermally stratified plasmas are unstable when the 
temperature decreases outwards, rather than when the en- 
tropy decreases outwards (the usual Schwarzschild criterion). 
This has been called t he magneto thermal instability (MTI). 
iParrish & Stond J2005I) show that in nonrotating atmospheres 
the MTI leads to magnetic field amplification and efficient 
heat transport. In future global simulations of RIAFs, it will 
be interesting to explore the combined dynamics of the MTI, 
the MRI, and angular momentum transport via anisotropic 
pressure stresses. 
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APPENDIX 
APPENDIX 

Grid and Variables 

Figure |A9] shows the location of variables on the grid. Scalars and diagonal components of second rank tensors (p, p\\, and 
p±) are zone centered. Vectors, representing fluxes out of the box, are located at the cell faces (V, B, and q||,±). The inductive 
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electric field (E) is located at cell edges such that the contribution of each edge in calculating ^ E • dl over the whole box cancels, 
and V • B = is satisfied to machine precision. The off diagonal part of the pressure tensor in Cartesian coordinates is related 
to n = bb(/?|| -p±)- This is a symmetric tensor whose components Pxy, Pxz, and Pyz are located such that the finite difference 
formulae for the evolution of velocities due to off diagonal components of stress are given by 

Vx,j,r' = Vx,j/-^(Pxylj^,,-Pxyij,)-^Pxz!!j^,,, -Px^j.t), (Al) 

Vzij.r' = Vzij/-^(Pxzl,j^,-Pxzlj^,)-pPyzlj^,,,-Py£!,j^,). (A3) 

Determination of 5t: stability, and positivity 

A time explicit algorithm must limit the time step in order to satisfy the Courant-Friedrichs-Levy (CFL) stability condition. 
Physically, St must be smaller than the time it takes any signal (via fluid or wave motion) to cross one grid zone. There is also a 
limit imposed on St for numerical stability of the diffusive steps. Additionally, since there are quantities which must be positive 
definite (p, p±), we also require St to satisfy positivity. We adopt the following procedure to choose St: 

_ mm{Sx,Sy,Sz} 

{\v\+\v,\ + m + \nL,\y ^^^^ 

min{Sx\Sy\Sz^ ^ (A5) 

2K|| 

min{Sx\Sy^,Sz^} 

Sti_ = , (A6) 

2k± 



where Va= B/ \/47r is the Alfven speed, and = max{ ^Jyp^Jp, ^/ip^Jp} is the maximum sound speed, taking the anisotropy 
into account. Stgdv, St\\, and and 5t±^ correspond to limits on the time step for stability to advection, and parallel and perpendicular 
heat conduction, respectively. 
The source steps for p\\ and p±^ are given by 

-L_JI = (^-V-q||-27?||b-VV-b+2^^V-^)"=Al, (A7) 

^-^-^^ = (-V-qxT-piV-V + pib-VV-b-^iV-/?)" =A2, (A8) 

where q^x = -^^x V|| T± denotes the temperature gradient part of qj^. For positivity of p"^^ and p'^' we require that the following 
conditions are satisfied: whenever A 1 andA2 are negative, St pos = min{-p||/Al,-p'^/A2}; if Al > 0, A2 < 0, ihsn5tpos = -p"^/A2\ 
if A 1 < 0, A2 > 0, then St pos = 'j /A 1 . Thus, our final constraint on the timestep St is given by 

(5f = Co X min|l/[max{(5fJ,, + ^f||^ + ^f2^}]'''^,min{(5fpoj}| (A9) 
where the max and min are taken over all zones in the box and Co is a safety factor (Courant Number) which we take to be 0.5. 

Implementation of the Pressure Anisotropy "hard wall" 

If the pressure anisotropy is larger than the constraints given in §2 by equations ( I32> - (I34> . then microinstabilities will turn on 
that will enhance the pitch-angle scattering rate and quickly reduce the pressure anisotropy t o near margina l stabil ity. Because 
this is a numerically stiff problem, we use an implicit approach, following the treatment of iBirn & Hessd J2001I) . Whenever 
equation (13 2> is violated, we use the following prescription for pitch angle scattering: 

pr=p\-\-p^^'-V-p"^'-i]^ (Aio) 

pT=P\ + ^^p5t [^-^-PT-'^j . (All) 

where Vp is a very large (^ 1 / ) rate at which marginal stability is approached. This implicit implementation (which can be 
solved by inverting a 2 x 2 matrix) with large Vp ensures that each time step the pressure anistropy will drop to be very near 
marginal stability for the firehose instability to break /i invariance. Given this pitch an gle sc atteri ng, th e collisionality parameter 
Veff in the thermal conductivity (eqs I42I - II44I ') is obtained by comparing equations (|A1Q> and JAl 1> with equations (I38> and 
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The effective pitch angle scattering rate v^ff is independent of Vp (and much smaller than Vp) in the limit of large Vp, and is by 
definition just large enough to balance other terms in equations ( I38I39> that are trying to increase the pressure anisotropy beyond 
marginal stability. 

The prescriptions for pitch angle scattering due to mirror modes and ion cyclotron waves are similar. For mirror modes we use 



pX'=p\-\.p5t[pX'-pT^2^-^„,^ 



«+i 



p'r=p\ + ^vp5t U«--;y-'+2C^ 



to limit the pressure anisotropy = 3.5 for our fiducial run ZlA) and Veff is given by 



'eff- 



For ion cyclotron pitch angle scattering we use 



: max < v„ 



(A13) 
(A14) 



(A15) 



pT=p"\\-1,^p^' 



<^-pT+s- 



(A16) 



and Veff is given by 



^n+\ „n+\ 



p"l' + S 



■■ max < 



\P\\ P± 



(A17) 



(A18) 



Implementation of the Advective part ofV ■ qj^ 



The flux of p±, q± = q±b, is given by 



Pl- 



iP\\-P^) 



B - VB 



B2 



P^ 



-KJ_V|| ( — ) +V,nagPL 



(A19) 



where the quantity in square brackets can be thought of as an advection speed due to parallel magnetic gradients. Because of 
this term, qj^ is not a purely diffusive operator, but also has an advective part characterized by the velocity V„,ag- If one treats the 
advective part via a simple central difference method, it does not preserve monotonicity. Instead, to treat the advective part of q^ 
properly, we include the advective part in the transport step. After including the advective heat flux in the transport step, it takes 
the form 

^ + V- [(V + y„,„,b)/.J =0. (A20) 

Thus, for updating p^ in the transport step we calculate fluxes on the cell faces using V + V,„flgb instead of just V. The transport 
step is then di rectionally split in the thre e directions. The procedure for monotonicity preserving schemes for calculating fluxes 
is described in lStone & Normal] ^1992a^ . 



Error Analysis 

The standard errors in the time averages reported in Table IA2I and in Figure IA7I are estimated by taking into account the 
finite correlation time for the physical quantities in the simulation, using techniques recommended by Nevins (2005). That 
is, the standard error for the time average {x) = J dtx{t)/T of a signal x{t) is given by cr^^^ = y/'Var(x)/Neff, where Var(;ic) = 
J dt{x(t)- {x))^/T is the variance of x, Neff = r/(2T,„,) is the effective number of independent measurements, T =15 orbits is the 
averaging time, and r,„, is an estimate of the integrated autocorrelation time. There are significant subtleties in determining th e 
integrated autocorrelation time from data. To deal with this, we use a windowing technique as recommended bv lNevinsI J2005I) . 
using Tin, = Jq dTC(T)W(T/T„), whcrc C(r) is the 2-time correlation function from the data, W(t/t„) is a smooth window 
function that cuts off the integral at r '--^ t„, and r„, ~ \^TTi„, (this gives results insensitive to the choice of window width for 
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Tint ^ T). IWinters et alJ J2003h found from comparing 3 realizations of shearing box MRl simulations that the magnetic stress 
had a variation of approximately ±6.5% after averaging over 85 orbits. The simulations we show here were averaged over 15 
orbits, so extrapolating from Winters et al. ( 2003) one might expect the uncertainties to be larger by a factor of w -^85/15 w 2.4. 
This is consistent with the typical error bars we report in Table lA2l and Figure lATl 
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Fig. Al. — Normalized growth rate (j/Cl) of the MRI versus normalized pressure anisotropy, (j>±-p±)/p\\ for /9 = 100, ^2^42/^^ = ^15/16, and two different 
kRS. Note that even a small anisotropy can stabilize the fastest growing MRI mode. The growth at large pressure anisotropy for kR=/Ois due to the mirror mode. 




orbits 

Fig. A2. — Time evolution of volume-averaged magnetic energy for the fiducial run Z/4. Time is given in number of orbits. There is a small decrease in the 
magnetic energy at fs 2 orbits when the pressure anisotropy is sufficient to stabiUze the fastest growing mode. However, small-scale kinetic instabUities limit the 
magnitude of the pressure anisotropy, allowing the magnetic field to continue to ampUfy. As in MHD, there is a channel phase which breaks down into turbulence 
at f» 4 orbits. 
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Fig. A3. — Time evolution of volume-averaged magnetic and kinetic energies, Maxwell, Reynolds, and anisotropic stress, and pressure (py : solid line, p±: 
dashed line) for the fiducial model Z/4. Time is given in orbits and all quantities are normalized to the initial pressure po- SVy = Vy+(i/2)Clxitad Ap = -p±). 
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Fig. A4. — Time evolution of volume-averaged pressure anisotropy (47r(p|[ -p±)/B^: soUd Une) for model Z14. Also plotted are the "hard wall" limits on 
the pressure anisotropy due to the ion cyclotron (dot dashed line) and mirror instabilities (dashed line). Ion cyclotron scattering is generally more efficient in the 
steady state. The limits on pressure anisotropy are applied at each grid point while this figure is based on volume averaged quantities. 
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Fig. A5. — Time evolution of volume-averaged pressure anisotropy (47r(p|| -px)/B^: solid line) for model ZZ8. Also plotted are the "hard wall" limits on the 
pressure anisotropy due to the ion cyclotron (dot dashed line) and mirror instabilities (dashed Une), although the ion cyclotron scattering limit is not applied in 
this simulation. The vohne averaged pressure anisotropy saturates at smaller anisotropy than the mirror threshold at ^ = 3.5, which is the only limit on pressure 
anisotropy used. 
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Fig. A6. — Time evolution of volume-averaged magnetic energy (dashed hne: Bf/87rpo. solid hne: B^/Snpo, dot dashed Une: By/Sirpo), total stress (Wxy/po) 
in units of 10"^, and pressure anisotropy for model Zll. Time is given in orbits and all quantities are normalized to the initial pressure pq. SVy = y3, + (3/2)n;c 
and Ap = (j>\\-p±). In this calculation there is no heat conduction and no isotropization of the pressiffe tensor. All resolved MRI modes are thus stabilized by 
pressure anisotropy and the 'saturated' state is linear anisotropic Alfv^n waves with no net angular momentum transport. 
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Fig. A7. — Space and time avera ge of the total s tress ((Wij/ Po)) versus l/(ki^Sz) for different runs. Error bars shown are based on estimates of the correlation 
time of the fluctuations described iri lNevinJ 
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Fig. A8. — Maxwell (qm: squares) and anisotr opic stress (a^: triangles) plotted against the collision frequency nomialized to rotation frequency (u/Q,). 
Transition to MHD occurs for > 30 (see Table lA3l . 



18 



Sharma et al. 



g\\x, g±x)i,j,k 
• 



P ** 




(^z> B^, q\\^, q±z)ijXi 



p 



(d, p\\, p±)ij,k 



-(Vy, By, q\\y, q±y\j,k 



X 
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Simulation parameters 
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Note. — Vertical field simulation with initial (i = 400. Z indicates that all simulations start with a vertical field, 7', 7?' indicate low and high 
resolution runs respectively. Z/4 is the fiducial run. Z/1, Z/il are the runs in CGL limit. ZM/ andZMh are the MHD mns. 

'^Wavenumber parameter used in Landau closure for parallel heat conduction (eqs. |13| & 1141 '). 

'^Imposed limit on pressure anisotropy for pitch angle scattering due to mirror instability (eq. 1331 ). Excluding Z/1, Zhl, and Z/8 all of these 
calculations also use a pressure anisotropy limit due to the ion cyclotron instability (eq. 1341 '). 
^ {{)) denotes a time and space average taken from 5 to 20 orbits. 

*'Ap = (p\\ -p±) 

^ T hese c ases run for only i=s 4 orbits at which point the time step becomes very small because regions of large pressure anisotropy are created 
(see l4l4l . 
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TABLE A2 
Statistics for Model Z/4. 



Quantity / 



((/)> 



{{Sf)) 



■2U1/2 



min(/) max(/) 



Sttpo 
!^ 

p6V^ 
PR 

Po 

(P\\~P±) B,Sy 

4ir(p|i-;)^) 



(B2/2) 

47r(p||-;j^) B^Bj. 

]? (s2/2) 
Wit 
(S-/87r) 

BPO 



0.083 


0.092 


0.016 


0.021 


0.662 


0.276 


0.318 


0.048 


0.036 


1.987 




017 


0025 


0032 


144 


0.102 


0.094 


0.014 


0.0184 


0.63 


0.125 


0.079 


0.0127 


0.715 


0.0264 


0.037 


0.034 


0.0032 


0.008 


0.348 


0.229 


0.277 


0.0434 


0.037 


1.856 


0.097 


0.113 


0.0147 


-0.072 


0.6211 


0.198 


0.129 


0.0178 


0.017 


0.654 


-1.366 


0.51 


0.098 


-2.632 


-0.083 


0.5895 


0.1043 


0.0067 


0.3744 


0.8611 


0.3323 


0.2725 


0.017 


-0.5307 


1.2704 


0.7356 


0.3718 


0.0714 


0.032 


1.807 


1.6574 


0.6598 


0.084 


0.4364 


3.7159 


0.5357 


0.3975 


0.024 


-0.9105 


2.084 


1.2287 


0.5504 


0.119 


0.0854 


2.7243 


0.99935 


2.3 X 10-5 


1.1 X lO""' 


0.9993 


0.9994 


3.557 


1.665 


_a 


1.1178 


7.929 


.144x 10' 


3.49 X 103 




585.4 


1.993 X 10"* 



^We calculate the enor using the autocorrelation time only for quantiti es that saturate to a steady state after 
5 orbits. Estimate for con'elation time r^,,, is based on the discussion in lNevini 120051) . p± and /7|| show a 
secular growth with time, so this way of expressing them as an average and an error is not applicable. 



TABLE A3 

Simulations with an Explicit Collision Term 





{{AttAp/B^)) 




(((^» 


^^^^ PO 









-1.41 


0.18 


0.082 


0.196 


1.09 


1 


1 


-1.47 


0.152 


0.072 


0.173 


1.14 


0.88 


3 


-1.43 


0.178 


0.08 


0.181 


1.02 


0.92 


10 


-1.35 


0.165 


0.071 


0.159 


0.96 


0.81 


20 


-1.24 


0.174 


0.070 


0.136 


0.78 


0.69 


30 


-1.01 


0.213 


0.070 


0.113 


0.53 


0.58 


40 


-0.87 


0.239 


0.070 


0.095 


0.4 


0.48 


100 


-0.43 


0.223 


0.06 


0.032 


0.14 


0.16 



"Ap = (p|| -p±) 



